Single‐cell analysis revealed a potential role of T‐cell exhaustion in colorectal cancer with liver metastasis

Abstract Liver metastasis (LM) is an important factor leading to colorectal cancer (CRC) mortality. However, the effect of T‐cell exhaustion on LM in CRC is unclear. Single‐cell sequencing data derived from the Gene Expression Omnibus database. Data were normalized using the Seurat package and subsequently clustered and annotated into different cell clusters. The differentiation trajectories of epithelial cells and T cells were characterized based on pseudo‐time analysis. Single‐sample gene set enrichment analysis (ssGSEA) was used to calculate enrichment scores for different cell clusters and to identify enriched biological pathways. Finally, cell communication analysis was performed. Nine cell subpopulations were identified from CRC samples with LM. The proportion of T cells increased in LM. T cells can be subdivided into NK/T cells, regulatory T cells (Treg) and exhausted T cells (Tex). In LM, cell adhesion and proliferation activity of Tex were promoted. Epithelial cells can be categorized into six subpopulations. The transformation of primary CRC into LM involved two evolutionary branches of Tex cells. Epithelial cells two were at the beginning of the trajectory in CRC but at the end of the trajectory in CRC with LM. The receptor ligands CEACAM5 and ADGRE5‐CD55 played critical roles in the interactions between Tex and Treg cell‐epithelial cell, which may promote the epithelial‐mesenchymal transition process in CRC. Tex cells are able to promote the process of LM in CRC, which in turn promotes tumour development. This provides a new perspective on the treatment and diagnosis of CRC.

with LM originating from the right and left sides of the colon manifests a survival discrepancy, despite the fact that left-sided colon cancer normally has a higher prevalence of liver metastases. 7For CRC liver metastases, if detected in time, liver resection is the most effective treatment strategy, with a cure rate of up to 20% and more than 50% of the patients could survive for at least 5 years. 9CRC liver metastases have a relatively specific biology, which means that it is crucial to explore the biological mechanisms involved.
Conventionally, it is believed that tumours originating in epithelial tissues are susceptible to distal metastasis via lymph node metastasis. 10In mouse models, lymph node metastasis appears to assist distal metastasis of tumour cells. 11Accumulating evidence indicate that the interactions between epithelial cells and T cells play a critical role during LM in CRC.During tumour progression, tumour-associated macrophages (TAMs) promote epithelial-mesenchymal transition (EMT) via the TGFα signalling pathway, driving tumour cells to acquire mesenchymal properties and migrate to new metastatic sites. 12TGFα promotes angiogenesis, tumour growth and metastasis in CRC, especially in LM. 13 These findings suggest that epithelial cells can regulate cancer cell LM through EMT transformation as well as the interaction with TAMs during the development of LM in CRC.In addition, tumour tissues also have abundant activated CD4 + T lymphocytes and CD8 + T lymphocytes during LM in CRC.In the presence of activated T cells, the killing effect of tumour cells becomes 'selective' and tumour cells could escape the killing action produced by T cells. 14Applying advanced next-generation sequencing techniques, single-cell RNA sequencing (scRNA-Seq) enables the exploration of gene expression patterns in individual cells, unveiling the diversity within cell populations. 15ScRNA-seq has been used to characterize CRC heterogeneity to develop a risk model for cancer based on mitochondrial autophagymediated profiling. 16Also, Xu et al. used scRNA-seq and observed a marked decrease in activated B cells in CRC with LM. 17

| Raw data
The GSE231559 scRNA-seq dataset was downloaded from the Gene Expression Omnibus website (https:// www.ncbi.nlm.nih.gov/ geo/ ).After removing mitochondrial, ribosomal and haemoglobin genes, the downstream data were read using the Read10X function of the Seurat package, 18 and the SCTransform function was used for data normalization.Batch effect was removed by the harmony package 19 (max.iter.harmony= 20, λ = 0.5).Next, based on the first 20 principal components, Uniform Manifold Approximation and Projection (UMAP) was performed for dimensionality reduction, and finally the cells were clustered into groups using the FindNeighbors and FindClusters functions at resolution = 0.1.Cell subpopulations were annotated using the marker genes provided by the CellMarker database (http:// bio-bigda ta.hrbmu.edu.cn/ CellM arker or http:// 117.50. 127. 228/ CellM arker/ ). 20

| Pseudo-time analysis
To characterize the mechanism of LM development in CRC, pseudo-time trajectory analysis was conducted using Monocle. 21fferentially expressed gene set (|log2FC| >0.25, min.pct= 0.25) between different groups was compared using the FindMarkers function to construct differentiation trajectories.We used the new-CellDataSet to create the cds object and filter low-quality cells, followed by using reduceDimension function to reduce data size by running the DDRTree algorithm.The cells were then ordered by orderCells function and the trajectory plot was produced by the plot_cell_trajectory function.The plot_genes_in_pseudotime function showed a scatter plot of the correlation between the expression of genes of interest and pseudotime.The branch involving more normal cells was taken as the beginning and end of the trajectory.

| Single-sample enrichment analysis
The Hallmark dataset was downloaded from The Molecular Signatures Database (MSigDB, https:// www.gsea-msigdb.org/ gsea/ msigdb).Based on the function of GSVA package, 22 the pathways related to cancer progression were extracted and single-sample enrichment analysis (ssGSEA) was performed.

| Cell communication analysis
To explore the types of interactions between cells in primary CRC and CRC with LM, cell communication analysis was conducted using the CellChat package. 23The probabilities of different cell types to produce signalling and cell-cell contact were inferred based on the CellChatDB database and overexpressed ligands or receptors in these cell types were identified.

| Statistical analysis
We used the Wilcox test to compare the differences in continuous variables between the two groups.All calculations were performed in the R language (version 4.3.1).A p < 0.05 was considered statistically significant.

| Cell types in primary CRC samples and CRC samples with LM
To investigate cell populations in primary CRC and CRC samples with LM and relevant molecular features, these two types of samples were included in our scRNA-seq data analysis.After clustering analysis, UMAP dimensionality reduction analysis and marker gene annotation were performed using the Seurat package.A total of 112,586 cells were grouped into nine cell subpopulations, namely, T cells (marked with GZMA, CD3G, CD3E and CD3D), macrophages (marked with FCGR3A, CSF1R and CD163), epithelial cells (marked with KRT8, KRT19, KRT18 and EPCAM), B cells (marked with MS4A1 and CD79B), plasma cells (marked with JSRP1, POU2AF1 and TNFRSF17), myofibroblast (marked with TAGLN, THY1, FAP, COL3A1 and COL1A1), mast cells (marked with TPSB2, HPGDS, CPA3 and SLC18A2), dendritic cells (marked with PTCRA, LILRA4, SCT and CLEC4C), endothelial cells (marked with CDH5, PECAM1, VWF, ADGRL2 and FLT1) (Figure 1A-C).The proportions of T cells and epithelial cells varied considerably in the two types of CRC samples.The proportion of T cells increased in CRC with LM (CRC: 51.8%, LM: 73.36) and that of epithelial cells decreased in CRC with LM (CRC: 16.2%, LM: 6.24%) (Figure 1D), indicating that T cells and epithelial cells played a crucial role in the process of LM in CRC.

| T-cell subpopulations in primary CRC and LM tumour samples
Subpopulation T cells was categorized according to the difference in the proportion of T cells in primary CRC and LM samples.Based on the expression annotation of marker genes, T cells were subdivided into natural killer (NK)/T cells, regulatory T cells (Treg) and exhausted T cells (Tex) (Figure 2A).Specifically, we found that the proportion of NK/T cells was higher in LM samples, while that of Tregs and Tex was higher in primary CRC samples (Figure 2B).Cytotoxic factors, including CST7, GZMA, GZMB, IFNG and NKG7, were upregulated in Tex_LM, Tex_CRC, NK/T cells_LM and NK/T cells_CRC (Figure 2C).2F).The expression levels of cell adhesion and apoptosis-related genes, cell proliferation-related genes and cell migration-related genes were also evaluated in T-cell subpopulations in primary CRC and LM tumour samples.It was found that the expression of cell adhesion and apoptosis-related genes, cell migration-related genes was upregulated in NK/T cells in LM samples (Figure 2G).The expression of genes related to cell adhesion, apoptosis and proliferation was upregulated in Treg in primary CRC samples (Figure 2H), while the expression of genes related to cell adhesion, proliferation (CORO1A) and migration (PFN1) was upregulated in Tex in LM (Figure 2I).These results suggested a possible involvement of these genes in the development of LM in CRC patients.

| Epithelial cell subpopulations in primary CRC and LM tumour samples
Epithelial cells could be categorized into six subpopulations, namely, epithelial cells 1, epithelial cells 2, epithelial cells 3, epithelial cells 4, epithelial cells 5 and epithelial cells 6 (Figure 3A).Epithelial cells 1 and 2 had higher proportions, with LM samples having more epithelial cells 1 and primary CRC samples having more epithelial cells 2 (Figure 3B).Analysis on the cellular functions showed that epithelial cells 1 were mainly associated with negative regulation of protein activity, transport, for example, peptide transport, negative regulation of proteolysis, negative regulation of hydrolase activity, negative regulation of peptidase activity (Figure 3C), while epithelial cells 2 were mainly associated with epithelial cell migration, epithelial mesenchymal transition, for example, mesenchymal cell proliferation, epithelial cell migration, epithelial to mesenchymal transition, negative regulation of cell adhesion (Figure 3D).These results were confirmed by gene expression patterns in epithelial cells 1 and epithelial cells 2 because protein metabolic activity-related genes, and transporterrelated genes were expressed at higher levels in LM (Figure 3E), while epithelial cell migration-related genes, and epithelial mesenchymal transition-related genes were expressed at higher levels in primary CRC (Figure 3F).To conclude, epithelial cells may play an important role in the process of EMT to promote metastasis in CRC.

| Pseudo-time analysis for Tex cells
Pseudo-time analysis was performed to characterize the transformation of Tex cells and Tregs during LM, and the results revealed two evolutionary branches (cell fate 1 and cell fate 2) of Tex cells during the transformation of primary CRC to LM (Figure 4A).As cell fate 2 was more abundant in LM samples, therefore we further focused on the changes in the expression of tumour development-related genes in cell fate 2. The results showed that the expression of CACYBP and HSPB1 in cell fate 2 gradually increased with pseudo-time progression (Figure 4B-D).Correspondingly, the expression levels of CACYBP and HSPB1 were also higher in LM samples than those in primary CRC samples (Figure 4E).

| Pseudo-time analysis for epithelial cells 2
During tumour metastasis to the liver in CRC, epithelial cells 2 were mainly associated with epithelial cell migration and epithelial mesenchymal transition.Therefore, to investigate the alteration of epithelial cells 2 during LM, pseudo-time analysis was performed.
The results showed that epithelial cells 2 in CRC were at the beginning of the trajectory and epithelial cells 2 in LM were at the end of the trajectory.During the differentiation process, there were two differentiation trajectories (FATE 1 and FATE 2), and eventually part of the epithelial cells 2 resulted in the phenomenon of LM in CRC (Figure 5A).The genes in the cells at the end of each branch were analysed by GO_biological process, and we found significantly enriched epithelial cell migration, epithelial cell proliferation, response to oxidative stress and regulation of T-cell activation (Figure 5B).We also examined the differential expression of genes related to cell adhesion, migration, and mesenchymal transition in primary CRC and LM.The results demonstrated that the expression levels of CD24, CEACAM1, EPCAM, MARVELD3, RHOB, SDCBP and SNAI1 were higher in FATE 2, while CD24, CEACAM1 and MARVELD3 were expressed at higher levels in LM (Figure 5C).The expression trend of genes related to cell proliferation, and differentiation was also consistent, as we found that expression levels of AREG, FOXA2, FOXO3, MYC, OVOL1, RGCC, S100P and SOX4 were higher in FATE 2, while SOX4, S100P and AREG were expressed at higher levels in LM (Figure 5D).

| Functional and cellular communication information on the involvement of epithelial cells 2 and Tex cell fate 2
Analysis on the tumour metastasis-related pathway activities of epithelial cells 2 and Tex cell fate 2 in patients with M0 and M1 stages showed that Tex cell fate 2 and Treg cell fate 2 might play important roles in metastatic tumours (Figure 6A).In epithelial cells 2, biological processes such as cell adhesion, mesenchymal transition, proliferation and migration played more important roles in metastatic tumours (Figure 6B).

Cell communication was analysed to further explore the interaction between Tex and Treg cells and epithelial cells 2 in promoting
EMT.The results showed that in cell-cell contact, CD8A-CEACAM5 and ADGRE5-CD55 receptor ligand pairs had a greater possibility of interaction between Tex cells and epithelial cells 2 (Figure 7A), and CD8A, CEACAM5, ADGRE5 and CD55 exhibited a greater possibility of interaction between epithelial cells 2 and Tex cells (Figure 7C).
In secreted signalling, the GZMA-F2RL1 receptor-ligand pair showed a greater possibility of interaction in Tex-epithelial cells 2 (Figure 7B) and was high-expressed inepithelial cells 2 and Tex (Figure 7D).

| DISCUSS ION
In this study, we analysed scRNA-seq data from primary CRC samples and CRC samples with LM.By constructing a single-cell atlas, it was found that the proportion of T cells in the LM increased and that of epithelial cells decreased.Study shows that that liver sinusoidal endothelial cells capture cancer cells, while macrophages phagocytose cancer cells and release tumour-killing cytokines. 24During this process, antigen-presenting cells present antigens to T cells for conversion into effector T cells, and this process could be enhanced by CD4 + T cells.However, when cancer cells evade the immune system, effector T cells under the effect of immune checkpoints would become dysfunctional. 24In cases of preoperative chemotherapy for primary CRC, the abundance of T cells reduces significantly. 25EMT is a process during which epithelial cells lose the characteristic features and turn to malignant tumour cells with metastatic tendency, which could explain the decrease in the number of epithelial cells.Studies have increasingly confirm that Tex cell formation undergoes metabolic deficiencies and is accompanied by overall alterations in signalling molecular pathways and epigenetics that could result in the suppression of effective immune responses and ineffectiveness of immune checkpoint therapies, partly contributing to tumour progression. 27,28There were three molecular subtypes of T cells (Treg, Tex and NK cells), with Treg and Tex as the predominate subtypes.In the tumour mesenchyme, chemokines produced by malignant cells and some other cells attract Tregs from the circulation to reach the tumour site.Pre-existing Tregs can be clonally amplified in response to specific antigen activation in the presence of TGFβ and IL-10, two abundant cytokines in the TME.These cytokines together with suboptimal antigen presentation promote the conversion of conventional T cells into Tregs with suppressive functions. 29Environmental factors also play an important role, and CRC often exhibits persistent antigenic stimulation and chronic inflammation in the TME.Chronic antigen exposure could increase excessive and sustain T-cell activation, eventually leading to T-cell depletion and the formation of Tregs. 30LAR and TNFRSF1B in Treg were high-expressed in primary CRC.Treg cells have a higher apoptosis rate compared to conventional T cells, which is correlated with low c-FLIP expression.Treg. 32These studies suggest that high expression of CFLAR and TNFRSF1B in Treg cells plays an important role in the regulation of the TME and progression of CRC, especially in regulating immune responses and promoting immune escape.There were six cell subpopulations in epithelial cells, and epithelial cells 1 and epithelial cells 2 were the dominant cell types.
Epithelial cells 1 was largely involved in the negative regulation of protein activity, and transport, while epithelial cells 2 were more associated with epithelial cell migration, epithelial mesenchymal transition.In cancers, the regulation of protein activity and transport in epithelial cells is closely related to cancer progression.PEPT1 and PEPT2 are important transporter proteins for the uptake of peptides and amino acids in intestinal renal epithelial cells, which regulate the uptake and metabolism of oxidative drugs and influence the therapeutic outcome of cancer. 33,34Epithelial cells may undergo significant structural changes in the complex cytokine and chemokine environment within the TME, which may lead to loss of cellular function. 35,36Epithelial cells 2 were mainly associated with epithelial cell migration, epithelial mesenchymal transition and were proportionally higher in primary CRC and had higher expression levels of BTG1, IER3, ZFP36L1 and DUSP1.Upregulated BTG1 expression can promote CRC cell proliferation. 37Blocking the IER3/ Nrf2 pathway significantly enhances the effect of antitumor drugs on bladder cancer. 38Significantly overexpressed ZFP36L1 in highly malignant tumour cells causes a high accumulation of proliferative and migratory properties in tumour cells. 39DUSP1 is significantly overexpressed in apatinib-resistant patients with gastric cancer, and downregulating DUSP1 may be a potential measure to attenuate apatinib resistance in gastric cancer. 40 Therefore, revealing the mechanism underlying LM occurrence based on scRNAseq technology might provide potential immunotherapeutic targets to facilitate the risk prediction for CRC patients.The current research analysed scRNA-seq data from primary CRC samples and CRC samples with LM to identify specific cell types and receptor/ligand factors involved in the interaction during LM in CRC.In addition, pseudo-time analysis was conducted to characterize the differentiation trajectory of epithelial cells and potential mechanisms of epithelial cells during LM.In conclusion, our study may provide a new direction for the treatment of CRC.
GO biological process analysis was performed based on the differentially expressed genes (DEGs).Specifically, NK/T cells were involved in the regulation of cell-cell adhesion, regulation of T-cell activation, regulation of inflammatory response and regulation of the apoptotic signalling pathway (Figure 2D); Treg was involved in regulation of T-cell activation, regulation of cell-cell adhesion and T-cell proliferation (Figure 2E); Tex was involved in regulation of cell-cell adhesion, regulation of T-cell activation and positive regulation of T-cell activation (Figure These results suggested that Tex and Treg cells played an important F I G U R E 1 Cell types in primary colorectal cancer (CRC) and CRC liver metastatic tumour samples.(A) Nine cell subpopulations in primary CRC and CRC liver metastatic tumour samples.(B) Bubble diagram showing the highly expressed marker genes in the nine cell subpopulations.(C) Violin plot demonstrating highly expressed marker genes in nine cell subpopulations.(D) The proportion of nine cell subpopulations in primary CRC and liver metastasis (LM) samples.F I G U R E 2 Landscape of T-cells subpopulations in primary colorectal cancer (CRC) and liver metastasis (LM) tumour samples.(A) Three subpopulations in T cells, NK/T cells, Treg, and Tex.(B) The proportion of NK/T cells, Treg, and Tex in primary CRC and LM tumour samples.(C) The expression levels of cytotoxic factor, exhausted factor, regulation factor, naïve factor, and co-stimulatory factor.(D) GO biological process involved in NK/T cells.(E) GO biological process involved in Treg.(F) GO biological process involved in Tex.(G) The expression levels of cell adhesion and apotosis-related genes and cell proliferation-related genes in NK/T cells.(H) The expression levels of cell adhesion and apotosis-related genes and cproliferation-related genes in Treg.(I) The expression levels of cell adhesion and profileraion-related genes and cell migration-related genes in Tex.role in epithelial cells during CRC development, indicating that Tex cells may be involved in tumour progression by altering the immune microenvironment in CRC.

F I G U R E 3
The subpopulations of epithelial cells in primary colorectal cancer (CRC) and liver metastasis (LM) tumour samples.(A) Six subpopulations in epithelial cells.(B) Proportion of the six subpopulations in primary CRC and LM tumour samples.(C) GO biological processes involved in epithelial cells 1. (D) GO biological process involved in epithelial cells 2. (E) The expression levels of protein metabolic activity-related genes and transporter-related genes.(F) The expression levels of epithelial cell igration-related genes, epithelial mesenchymal transition-related genes.Interactions between cells in the tumour microenvironment (TME) may increase the proportion of T cells and reduce the proportion of epithelial cells.The reduced proportion of these two types of cells may be responsible for the development of LM in CRC.

F I G U R E 4
Pseudo-time analysis for Tex cells.(A) Pseudo-time differentiation trajectory of Tex cells.(B) Heatmap of the expression of tumour development-related genes in the pseudo-time trajectory of Tex cells.(C) Expression heatmap of tumour development-related genes in cell fate 1 and cell fate 2 in the pseudo-time trajectory of Tex cells.(D) Expression changes of CACYBP and HSPB1 in cell fate 1 and cell fate 2 with pseudo-time progression.(E) Expression levels of CACYBP and HSPB1 in primary colorectal cancer (CRC) and liver metastasis (LM) samples.Specific knockdown of c-FLIP in Treg cells results in a fatal autoimmune disease that manifests as peripheral Treg cell deletion, effector cell activation, multi-organ immune cell infiltration and premature death. 31TNFR2 + Treg cells preferentially accumulate in the TME and express high levels of immunosuppressive molecules.In an in vitro model, the TNFα/TNFR2 pathway upregulates Foxp3 expression in CD4 + CD25 + T cells and latent TGFβ production in Treg and enhances the immunosuppressive function of F I G U R E 5 Pseudo-time analysis for epithelial cells 2. (A) Pseudo-time differentiation trajectory for epithelial cells 2. (B) Heatmap of tumour progression-associated GO terms in the pseudo-time trajectory for Epithelial cells 2. (C) Expression changes of genes related to cell adhesion, migration and mesenchymal transition with pseudo-time progression in FATE 1 and FATE 2, different samples.(D) The expression changes of cell proliferation, differentiation-related genes in the FATE 1 and FATE 2, different samples, with pseudo-time progression.
It should be noted that further studies are needed to elucidate the specific mechanisms underlying the high expression of epithelial cells 1 and epithelial cells 2 in CRC liver metastases.This study revealed multiple paired transcripts that mediated the communication between Tex and Treg cells and epithelial F I G U R E 6 Functions involved in epithelial cells 2 and Tex cell fate 2. (A) Single-sample gene set enrichment analysis (ssGSEA) enrichment scores for Tex cell fate 2 and Treg cell fate 2. (B) The ssGSEA enrichment score of Epithelial cells 2 in tumour metastasis-related signalling pathways.cells, such as ADGRE5-CD55 and CD8A-CEACAM5.CEACAM5 as a molecule closely related to cell adhesion 41 plays a key role in the therapeutic diagnosis and prognosis in a variety of epithelial tumours, including respiratory and gastrointestinal tumours. 42Gebauer et al. indicate that CEACAM5 overexpression can promote tumour growth through the EMT pathway in the localized regions of aggressive cancer. 43ADGRE5 (CD97) has been considered as the most promising target of adhesion G protein-coupled receptors in glioblastoma. 44ADGRE5 in retinal pigment epithelial cells also controls leukocyte activation and trafficking in uveal retinal inflammation. 45CD55, a complement regulatory protein, binds to multiple ligands that are constitutively expressed on monocytes or granulocytes and is rapidly upregulated during activation of T cells. 46It has been shown that CD55 regulates both innate and adaptive immune responses and mediates T-cell co-stimulation when binding to CD97. 47This revealed a potential interaction between Tex and Treg cells in the colon epithelium, which in turn promoted mesenchymal transformation and the development of cancer cells.However, there were some limitations in our study.First, the main source of data came from public databases, F I G U R E 7 Epithelial cells 2 and Tex cell fate 2 cells exchange information.(A) Receptor-ligand pairs in cell-cell contact, Tex and Treg in cell-cell contact with epithelial cells 2. (B) Secreted signalling, receptor-ligand pairs for cell-cell contact of Tex and Treg with epithelial cells 2. (C) The expression of receptor-ligand pairs for cell-cell contact of Tex and Treg with epithelial cells 2 in single-cell data.(D) The expression of receptor-ligand pairs of Tex and Treg for secreted signalling with epithelial cells 2 in single-cell data.